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Abstract  -  A  Monte-Carlo  simulation  technique  for  the  calculation  of  the 
partition  function  of  a  general  Gibbs  random  field  is  presented.  We  sho'™'  that 
the  partition  function  of  a  general  Gibbs  random  field  is  equivalent  to  an 
expectation.  This  observation  allows  us  to  develop  an  importance  sampling 
approach  for  estimating  this  expectation  by  using  Monte-Carlo  simulations. 
Two  different  methods  are  proposed  for  this  task.  We  show  that  the  resulting 
estimators  are  unbiased  and  consistent.  Computations  are  performed  itera¬ 
tively,  by  using  a  simple.  Metropolis-like,  Monte-Carlo  algorithm  with  remark¬ 
able  success,  as  it  is  demonstrated  by  our  simulations.  Our  work  concentrates 
on  binary,  second-order  Gibbs  random  fields  defined  on  a  rectangular  lattice. 
However,  the  proposed  methods  can  be  easily  extended  to  more  general  Gibbs 
random  fields  and,  therefore,  can  become  quite  useful  in  many  scientific  areas, 
such  as  biology,  statistical  mechanics  and  image  processing.  Furthermore,  a 
potential  contribution  of  our  technique  to  optimally  estimating  the  parameters 
of  a  general  Gibbs  random  field  from  a  given  realization  via  a  maximum- 
likelihood  approach  is  anticipated. 
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I.  INTRODUCTION 

Markov,  or,  equivalently,  Gibbs  random  fields  (GRF’s)  belong  to  a  well 
known  and  popular  class  of  parametric  random  field  models  [1],  [2].  They  are 
extensively  used  for  modeling  spatial  interaction  phenomena,  including  the 
phenomenon  of  phase  transition  which  occurs  when  some  macroscopic  proper¬ 
ties  of  a  physical  system  change  discontinuously  over  a  small  perturbation  of 
its  parameters.  Some  “classic”  application  areas  of  GRF’s  include  statistical 
mechanics  [3-5],  ecology  [6],  sociology  [7]  and  crystallography  [8],  [9],  Recently, 
GRF  models  have  found  wide  applicability  in  the  general  areas  of  image 
analysis  and  computer  vision,  where  they  have  become  an  attractive  tool  for 
providing  statistical  models  for  images  [10-13].  Since  they  capture  various 
image  characteristics  in  terms  of  a  few  parameters,  they  have  been  successfully 
applied  in  a  variety  of  image  processing  problems,  such  as  smoothing  and  seg¬ 
mentation  [14-18],  restoration  [19],  [20],  reconstruction  [21]  and  coding  [22],  as 
well  as  in  a  variety  of  computer  vision  tasks  [23],  [24].  However,  many  theoret¬ 
ical  and  computational  difficulties  prohibit  the  application  of  these  models  to  a 
wider  class  of  problems,  one  of  the  major  difficulties  being  the  computation  of 
the  partition  function.  No  exact  solutions  are  known  for  this  problem,  except 
for  very  simple  cases  [5],  which  are  not  adequate  for  most  applications  of 
interest.  As  an  alternative  to  exact  solutions,  approximation  methods  are  used 
to  get  a  good  estimate  of  this  quantity.  Some  of  the  better  known  approxima¬ 
tion  techniques  could  be  grouped  into  the  following  categories  [5]:  (a)  cell,  or 
cluster,  approximations  [25],  [26],  (b)  series  expansions  in  powers  of  an 
appropriate  variable  [27],  [28],  (c)  renormalization  group  techniques  [29],  [30]; 
and,  (d)  statistical  simulation  and  other  techniques  [5],  [31],  Most  of  these 
methods  are  either  unreliable,  especially  around  the  critical  point  where  a 
phase  transition  occurs,  or  require  considerable  faith  in  the  involved 
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assumptions  [5]. 

The  lack  of  a  closed  form  solution  for  the  partition  function  imposes  many 
restrictions.  For  example,  the  statistical  inference  of  GRF’s,  via  a  maximum- 
likelihood  approach,  is  an  open  problem,  since  solution  of  the  likelihood  equa¬ 
tions  is  impossible  without  knowing  the  exact  value  of  the  partition  function 
and,  possibly,  some  of  its  derivatives.  As  a  result  of  this,  alternative  methods 
have  been  developed  with  moderate  success  [6],  [11],  [17],  [20],  [32-36].  It  is, 
therefore,  clear  that  the  problem  of  calculating  the  partition  function  of  a  gen¬ 
eral  GRF,  in  a  computationally  amenable  way,  is  of  great  interest.  The  study  of 
this  problem  is  the  purpose  of  the  present  paper. 

The  approximation  technique  proposed  here  is  based  on  a  stochastic  simu¬ 
lation  approach  and  it  makes  use  of  a  mutually  compatible  Gibbs  random  field 
(or,  equivalently,  Markov  mesh  random  field)  and  its  relation  to  a  general  GRF 
[37],  [38].  Section  II  is  devoted  to  establishing  the  required  background  and 
notation.  In  Section  III  we  discuss  the  problem  of  computing  the  partition 
function  of  a  general  GRF  via  Monte-Carlo  simulations,  and  we  propose  two 
methods  to  achieve  this,  which  correspond  to  different  approximations  of  a  gen¬ 
eral  GRF  by  a  MC-GRF  random  field.  In  Section  IV  we  consider  the  algo¬ 
rithmic  implementation  of  these  methods  and  discuss  various  important  pro¬ 
perties  of  the  obtained  estimators.  Estimation  of  the  derivatives  of  the  parti¬ 
tion  function,  and  other  related  quantities,  is  very  important  and  is  also  neces¬ 
sary  in  order  to  carry  out  our  second  method.  This  is  discussed  in  Section  V, 
whereas,  in  Section  VI  various  supporting  simulation  experiments  are 
presented.  Finally,  in  Section  VII,  we  review  our  results  and  draw  our  conclu¬ 


sions. 
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II.  GIBBS  RANDOM  FIELDS 

Assume  that  we  have  a  collection  of  MxN  sites  which  form  a  two- 
dimensional  rectangular  lattice  AMN  =  {( i,j ):  1<  i  <M  ,  1<  j  <N  ).  A  discrete  ran¬ 
dom  variable  H(i,j)  is  assigned  at  each  site  of  the  lattice,  taking  values  from  a 
discrete  ensemble  EH  =  fo,  .  . . ,  },  which  contains  R>  2  distinct  values.  The 

resulting  random  field  [H]  =  (H(i ,  j):  1<  i  <M  ,  1<  j  <N  }  can  take  any  one  of 
the  Rmn  possible  realizations  [A]  =  ( Ay  :  1<  i  <M ,  l<j  <N  }  in  the  Cartesian 
product  Ehn  with  joint  probability  distribution1  Pr [H  =  h],  At  each  lattice  site 
sn,  n  =  1,2, ,  MN ,  a  set  N3n  of  neighbor  sites  is  assigned  such  that:  (a)  sn  4 
iVJit  ;  and,  (b)  if  sneN,m,  then  smeNSri.  We  shall  restrict  [Pf]  to  be  in  the  class  of 
Markov  random  fields  (MRF’s).  In  this  case  the  following  Markovian  property 
is  satisfied: 

Pr[H(sn)  =  h(sn)\H(s)  -  h(s),  s  e  AMN  -  {s„  )  ] 

=  Pr[H(sn)  =  h(sr.)\ms)  =  h(s),  s  eN,n]  >  0. 

According  to  the  Hammersley-Clifford  Theorem  [6],  [H]  is  equivalent  to  a  GRF 
whose  joint  probability  distribution  is  given  by  the  Gibbs  measure 

Pr[H  =  h]  =  -|exp{--|rt/(h)}  .  (1) 

In  eq.  (1),  Z  is  a  normalizing  constant  known  as  the  partition  function,  T  is  a 
positive  parameter,  known  as  the  temperature,  which  controls  the  degree  of 
peaking  in  the  Gibbs  measure,  whereas,  U(-)  is  the  energy  function  which 
depends  on  a  specific  realization  h  of  the  GRF  and  the  potentials  associated 
with  its  cliques  [6]. 

1  In  the  following,  boldface  characters  denote  vectors  which  correspond  to  a  lexicographic 
ordering  of  the  random  field. 
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Without  loss  of  generality,  we  consider  second-order  GRF’s  [6].  These  ran¬ 
dom  fields  is  a  “rich”  class  of  models  capable  of  describing  a  variety  of  spatial 
interaction  phenomena  in  a  satisfactory  way.  In  this  case,  the  neighborhood 
Nif,  at  site  (i,j)e. AMN,  will  be  given  by 

N&  =  {(i-1,  y- 1),  (i-1,  j),  (i-1,  j+ 1),  (i ,  j~  1).  ( i ,  j+ 1),  (i  +1, 7-1),  (i +1,  j),  ( i  +1, y'+l)}  , 
and  the  Gibbs  measure  (1)  by 


M  N 


Pr[H  =  h]  = 


rr  n  (z*y  >  7  >  ^«-i,  y-i»  h «,y-i)  > 


(2a) 


i=l  7=1 


where 


M  N 


state 8  h  1=1  7=1 


(2b) 


In  eq.  (2),  atJ (x ,  y ,  z,  to)  is  the  ZocaZ  transfer  function  (LTF)  of  the  GRF  [ff],  and 
the  summation  is  carried  over  all  RMN  states  [37].  The  LTF  has  to  be  modified 
at  the  boundary  sites  of  AMN  depending  on  the  type  of  boundary  conditions 
assumed  (free,  or  toroidal).  In  this  paper,  we  shall  assume  free  boundary  con¬ 
ditions2;  i.e.,  H(i,  j)  =  <)>!  in  [  (1,7)  ,  i  <  0  or  7  <  0  }.  The  case  of  a  first-order 
GRF  with  neighborhood 

=  {(i-1,  y  ),(i ,  7-I), (i ,  7+1). (*+l,  7)}  , 

is  a  special  case  of  eq.  (2),  with  ai7(Ay,  Zx,_1>7-,  Ai)7_i)  not  depending  on 

A  special  case  of  a  general  GRF  is  a  mutually  compatible  Gibbs  random 
field  (MC-GRF),  or,  equivalently,  a  Markov  mesh  random  field  [37].  The  proba¬ 
bility  structure  of  a  subset  [H]A  of  such  a  random  field  [if],  restricted  on  a 
finite  sublattice  A  of  AMN,  is  independent  of  the  size,  or  shape,  of  A.  The  LTF 
x ij(x,y,  z,  < a)  of  these  GRF’s  is  restricted  to  satisfy  the  following  relationship: 


2  Thi3  type  of  boundary  conditions  is  most  natural  in  many  practical  situations  [19]. 
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^  ty  (u  ,y ,  z ,  a)  =  kij  , 
u<=Eh 

for  every  triplet  (y  ,z  ,a>)eEl31,  where  fey  is  a  constant.  It  is  easy  to  show  that  the 
computation  of  the  partition  function  ZMN  in  this  case  is  trivial.  Indeed, 

M  S 

z>mn  =  x  nn^.Ai-u.WxAy-i) 

states  h  i=l  j= I 

MS  MS 

=  nn  x  ty (u , j_i, ■ 

»=l7=l  »€£„  i=lj=l 


In  this  case,  the  partition  function  is  expressed  as  the  product  of  regular  and 
local  partition  functions  (the  fey’s);  therefore,  no  phase  transition  is  associated 
with  these  models  [37]. 

Generating  a  realization  of  a  MC-GRF  is  also  an  easy  task,  because  this 
can  be  done  lexicographically  using  point  by  point  simulation,  as  the  following 
relation  implies: 


M  N 

Pr[H  =  h  ]  =  nn 


»=i v=i 


tjj  (by ,  hj-i,j  thj-ij-uhj'  j~i) 

X  Ty  (u  >  h- 1,  j ,  hi- 1,  j- 1,  fe; _  j- j) 

U€E„ 


M  S 


—  h-ij  I  fej-i,  j ]  , 

«=w=i 


where 


Rr  [  fey  I  hi-i  j,hi-\j-\,hi  j-i  ]  — 


Xy  (fey  ,  fej  — 1_  j  ,  fej  — 1,  _/ — 1,  fej  ,  J  — l) 

X  Xy  (u  ,fe,'_l(  y  ,  hi -l  j-l,  hij-l) 

U€% 


Therefore,  knowing  the  values  of  random  variables  H(i-l,j),  H(i-l,j-l)  and 
H(i,j- 1),  we  can  generate  the  value  of  H(i,j)  by  sampling  the  conditional  pro¬ 
bability  Pr[  fey  I fej_i,  j , fe;-i,  j~i, hi  ]  over  all  possible  values  in  EH. 

Since  a  MC-GRF  is  characterized  by  many  attractive  properties,  it  has 
been  pointed  out  that  it  might  be  a  good  idea  to  try  to  approximate  a  general 
GRF  by  a  MC-GRF  [38],  Two  approaches  have  been  adopted  here  for  the 
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solution  of  this  approximation  problem.  The  first  one  is  the  simplest  one,  how¬ 
ever  the  second  one  achieves  optimality  with  respect  to  an  “entropy”  distance 
between  the  probability  distribution  of  a  general  GRF  and  the  probability  dis¬ 
tribution  of  a  MC-GRF.  Both  will  be  used  to  develop  an  importance  sampling 
procedure  for  the  Monte-Carlo  calculation  of  the  partition  function  ZMN,  given 
by  eq.  (2b).  This  will  be  discussed  in  the  following  section. 


III.  MONTE-CARLO  CALCULATION  OF  THE  PARTITION  FUNCTION 


As  we  mentioned  before,  the  main  purp  se  of  our  work  is  to  calculate  the 
partition  function  ZMN  of  a  general  GRF  via  Monte-Carlo  simulations.  From  eq. 
(2b)  observe  that 


Z'MN  -  X  A-mn  (h)  I 

states  h 


(3a) 


where 


M  N 

A-mnO*)  = 

.=t  jm  1 


(3b) 


Since  the  previous  summation  is  carried  over  RMN  states,  it  is  not  possible  to 
compute  ZMn  directly,  even  for  moderate  size  lattices.  Assume  that  P<Wv(h)  is  a 
joint  probability  distribution  defined  on  the  space  of  all  RMN  realizations  h, 
such  that  PiWV(h)  >  0,  for  every  state  h  e  Ehn .  In  this  case 


Zmn  =  X 

states  h 


Pwv(h) 


Pmn  (h)  =  EP 


Amn( h) 
Pmn( h) 


(4) 


From  eq.  (4)  we  derive  the  following  estimate  ZMN  P  for  the  partition  function 
Zmn- 


%mn,p  -  £ 


MN,p(K)  =  lim 

X-»+~ 


_1_  y  Amn  (hk  ) 

K&PMNi  h*) 


(5a) 


where 
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provided  that 


Zmn,p 


KAMN(hk) 

A=i-fW(h*) 


Amn  (h) 

-  <  +oo 

PMN(  h) 


(5b) 


(5c) 


for  every  h  e  EffN.  In  eq.  (5),  h* ,  &  =  1,  2 ,  K,  are  i.i.d.  realizations  drawn 

with  probability  Pun  (h*),  and,  therefore,  2mn,p  is  a  Monte-Carlo  estimate  of 
the  partition  function  [39].  Inequality  (5c)  is  satisfied  in  the  case  of  finite  lat¬ 
tices.  When  MN  +°°,  this  inequality  may  be  violated  and  our  approach  has  to 
be  modified.  This  modification  will  be  discussed  shortly. 

Let  us  now  focus  on  the  problem  of  choosing  the  appropriate  joint  proba¬ 
bility  distribution  PMN(h).  Since  ZMNiP(K),  and,  therefore,  ZMNP,  are  unbiased 
estimators  of  ZMN  [33],  *,ve  shall  concentrate  on  finding  joint  probability  distri¬ 
butions  PMN(h)  which  result  in  a  small  error  variance  EP[  (  ZMN  P{K)  -  ZMN  )2  ], 
for  every  K.  Ultimately,  we  would  like  to  minimize  the  quantity 


states  h 


■A&y(h) 
Pmn(  h) 


K  -  1 


with  respect  to  PMN( h),  provided  that  we  can  draw  samples  h  from  the  proba¬ 
bility  distribution  PMN(h)  in  a  computationally  efficient  way.  Equivalently,  wre 
would  like  to  minimize  (with  respect  to  PMN(h))  the  Ali-Silvey  distance  [40] 


d/s (kmn ,Pmn)~  X 

states  h 


rc,w/(h) 

^W(h) 


2 

iWb) 


=  -JI~EPl2fiN'P(K)]-(K-l)  ,  (6) 

Zmn 


between  the  probability  measures  Pmn  (h)  and  rt^Ch).  Notice  that  P,w;v(h) 
should  be  known  explicitly,  since  at  each  step  of  the  Monte-Carlo  algorithm  the 
quantity  A.MN(h)/  PMN(h)  must  be  calculated. 
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It  is  straightforward  to  show  that  the  direct  minimization  of  eq.  (6),  with 
respect  to  PWJV(h),  results  in  PMN( h)  =  rcWiV(h)  =  AMN(h)/ZMN  ,  which  is  the  Gibbs 
measure.  Although  this  choice  results  in  zero  error  variance,  it  is  useless,  since 
the  computation  of  AMN(h)/  PMN(h)  requires  ZMN  to  be  known  explicitly.  We  can 
resolve  this  problem  by  considering  probability  distributions  PMN(h)  which 
correspond  to  a  MC-GRF,  therefore,  satisfying  the  following  relation: 

M  N 

Pmn (h)  =  nil  *ij(?iij 'hi-ij .hi-ij-uhij-J  >  0  ,  (7a) 

i=ij=i 


where 


£  z,  <o)  =  1  ,  for  every  (y,z,a)<zE (7b) 

ueEH 

We  shall  denote  this  class  of  distributions  by  D;  i.e., 

iWh)  e  D,  if  and  only  if  eq.  (7)  holds  . 

This  is  a  convenient  choice,  since  drawing  samples  from  PMN{ h)  (i.e.,  generating 
a  realization  of  a  MC-GRF)  is  an  easy  task,  as  it  was  mentioned  in  the  previ¬ 
ous  section.  Computing  the  ratio  AMN(h)/PMN(h)  is  also  straightforward.  Indeed, 
if 


Gjj  (hjj  i  hj~  l.  j » ^j-i,  j-u  h-i,  j-0 
lijihij  >  l,y'-l,  ^i,  j-0 


(8a) 


and 


M  N 

Qmn,p(  h)  =  nn^j.p  (hij ,  hi_it  j ,  A,_i  j_i,  hi  j_0  ,  (8b) 

then  (see  eqs.  (3),  (4),  (7)  and  (8)) 

%mn=  2  Qmn.p^O  Pmn(&)  ,  (9) 

states  h 


and 


Zmn.p 


=  lim  ZMN  P(K)  =  lim 

K  —*  +°*  K —* 


S  Q.w.v.fih*)  , 

A  *=i 


(10a) 
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wheie 

%MN  ,p(K)  =  QMN.pfak)  »  (10b) 

K  k=l 

provided  that 


Qmn.rQ1)  <  +°°  >  (11) 

for  every  h  e  E^N.  The  error  variance  in  this  case  will  be  given  by  (see  eq.  (10;) 

VarP[ZMN  P  (K)]  =  —  £  Qmn ,  j>  (h) Pmn (b )  -  Z^  .  (12) 

A  states  h 

A  simple  choice  for  the  LTF  of  a  MC-GRF  could  be  xjjd(x,y  ,z,co)  -UR,  for 
every  (x  ,y  ,z ,a>)eE£  and  every  (i,j)e  AMN.  In  this  case  P^dN(h)=  URMN  (i.e.,  the 
joint  probability  distribution  of  an  i.i.d.  random  field)  and  the  error  variance 
will  be  given  by 

Varpiidl2UN'Pud{K)]  =  ±\RMN\  £  A^(h)l  -  Z^J  .  (13) 

states  h 


In  most  practical  situations  one  does  not  expect  the  i.i.d.  choice  to  give  a 
good  estimate  of  ZMN  in  a  reasonable  time.  The  reason  for  this  is  that,  in  most 
cases  of  interest  (e.g.,  in  low  temperatures),  only  a  small  fraction  of  realizations 
h*  contribute  substantially  to  the  partition  function  sum,  while  the  contribu¬ 
tion  of  most  other  realizations  is  negligible.  Considering  samples  h*  with  equal 
probability  will  greatly  underestimate  ZMN,  since  most  of  the  time  h*  will  be  a 
sample  that  does  not  contribute  much  to  the  computation  of  ZMN .  To  overcome 
this  difficulty,  we  will  need  many  iterations,  which  could  result  in  a  K  being 
comparable  to  RMN .  Nevertheless,  this  sampling  scheme  becomes  attractive  in 
the  case  of  GRF’s  which  are  “close”  to  i.i.d.  random  fields  (i.e.,  in  high  tempera¬ 
tures;. 

If  we  now  choose  the  joint  probability  distribution  of  a  MC-GRF  which  is 
“as  close  as  possible”  to  the  Gibbs  measure,  then  this  choice  will  favor  the  most 
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probable  realizations  of  the  GRF  over  the  less  probable  ones.  This  will  result  in 
the  error  variance  given  by  eq.  (12)  which  will  be  much  less  than  the  error 
variance  given  by  eq.  (13).  It  is  now  natural  to  seek  the  optimal  joint  probabil¬ 
ity  distribution  Pj^Ch)  of  a  MC-GRF,  by  solving  the  following  constrained 
minimization  problem: 

Pff/v(h)  =  arg  {  min  dls(nMN,  PMN)}  .  (14) 

Pmn  e  D 


However,  the  elution  of  eq.  (14)  is  not  feasible  in  general.  Therefore,  we  shall 
derive  two  suboptimal  solutions  to  this  problem.  The  first  one  is  motivated  by 
the  following  theorem: 

Theorem  T.  For  every  joint  probability  distribution  P,^(h)  e  D,  there  will 
be  a  joint  probability  distribution  P  ff^(h)  e  D,  with 

tljXx  ,y,z,<o)  =  xjjXx  ,y ,  z ,  co)  ,  (15a) 

for  every  (i,j)e  AMN  -  (( M,N ))  and 


tfih(x,y  ,z  ,'*) 


<yMN(x,y  ,z  ,a) 

Z  aww(u,y  ,z,u))  ’ 

U€% 


(15b) 


for  every  (x,y,z, co) e  E,),  such  that  Varpa)[ZMN  f  ^  Varpa)[2,MN  p(V(K)]  . 

Proof:  For  every  (y ,  z ,  w)  e  Ef],  let  *jfih(x,y,z,  co)  be  the  LTF  which  minim¬ 
izes  the  following  constrained  minimization  problem: 


minimize  £ 


CMN(u,y,z,<- a) 
iMN(u,y,z,  co) 


(16a) 


such  that 


Z  twv(u,y,z,a>)  =  1  ,  (16b) 

“f£« 


and 


*mn(x  ,y  ,Z ,  co)  >  0 ,  for  every  x  e  EH  . 


(16c) 


page  12 


By  usir  g  a  Lagrange  multiplier  kyia  and  by  differentiating  eq.  (16a)  with 
respect  to  xMN{x,y  ,z,  co),  we  obtain 


3 

dxMN{x,y,z,(o) 


X 

ue£* 


CMN(U,y,Z><») 


+  \as  [  X  -  1  1 

u<bEh 


=  0  , 


or 


*mn(x  >y  >z  ,(&) 


Gmn(x  ,y  ,z  ,co) 


>  o. 


(17) 


From  eq.  (16b)  we  have  that 


X  ,y,2,oj) , 

U€Eh 


which  together  with  eq.  (17)  gives  eq.  (15b).  Notice  that,  the  only  non-zero 
entries  of  the  Hessian  of 


X 

ueEH 


Wu,y,z,co)  + 


X  *MNto,y,Z 

ueEH 


co)  -  1 


will  be  the  diagonal  ones,  which  are  positive.  Therefore,  the  LTF  x'fihix  ,y  ,z ,  co), 
given  by  eq.  (15b),  is  the  solution  of  problem  (16).  If  \'MN  =  AMN  -  {(M ,  N)),  we 
have  that  (see  also  eqs.  (7),  (8),  (12)  and  (15)) 


K  Varr: 


■ft 


P<»  MN  ,Pa> 


IK) 


+  Z2 


MN  ~ 


staUshonAMN  i=l  j=l  ^y >  ^i-1,  j  >  ^i-l,  j-ls  ,j-0 


stales  hon  A' 


MN 


nn 


Qij^hij  A-i.yA-i,  j-i,  hij-i) 


o’,  j ,hj 


GmN^U'  i  ^i-1,  j  >  y-l>  ,j- 1) 
=£w  Ai-1,  j  ’  hi-ij-1,  ,  y-i) 


2  X 

states  bon  A’^y 


r-r  j— r  ®ij(h-ij  ,/li-l.y,  ^i-l.y-l,  hij-i) 

(, . ! )1a'^  t,f  (Ay ,  Ai-!.  j ,  Ai-1.  y-iA . >-i) 


X 


Gmn^u  A-i,;  A-i,y-i>  hij-0 


ugEh  >  h-i-l,j  i  ^i-1,  y-i>  ,  y-i) 


y  Gij(hjj  ,  A|-1,  j  ,  A|-1  y-1,  /t,  y-i) 

states  hon  Aww  i=l  j=l  ^y  '(^y  >  ^>-1,  j  >  ^i-1,  j-lt  Hi ,  y-l) 


=  KVar 


pi  2) 


2WA,,p(2)A) 


+  .  (18) 
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The  proof  is  now  a  direct  consequence  of  eq.  (18). 

□ 

Given  the  LTF  Ogix  ,y  ,z  ,m)  of  a  GRF,  our  first  suboptimal  choice  for  the 
joint  probability  distribution  PMN(h)  will  be  a  joint  probability  distribution 
pMN(h),  given  by  eq.  (7),  with  LTF 


x'jix ,y  ,z  ,(£>) 


Ogix  ,y  ,z  ,<s>) 

Z  csgiu  ,y  ,z  ,(a) 

U€Eh 


(19) 


for  every  (x ,  y ,  z ,  <o)  e  E £  and  every  ( i,j)e  AMN.  In  this  case  (see  eq.  (8a)) 


QijP‘0c,y,z,oi)=  X  dij{u.,y  ,z,(a) .  (20) 

ueEH 

The  MC-GRF  with  LTF  given  by  eq.  (19)  corresponds  to  the  MC-GRF 
obtained  by  approximating  a  general  GRF  via  the  “Approach  B”  developed  in 
[38] .  This  is  a  sensible  choice,  because  the  computation  of  the  joint  probability 
distribution  7W(h)  is  feasible,  and  the  obtained  MC-GRF  contains  substantial 
information  about  the  original  GRF  [38],  much  more  than  the  i.i.d.  random 
field,  which  contains  absolutely  no  information  about  the  GRF  in  low  tempera¬ 
tures.  One  expects  that  the  samples  h,  drawn  from  this  joint  probability  distri¬ 
bution,  will  contribute  substantially  to  the  computation  of  the  partition  func¬ 
tion  sum,  resulting  in  a  considerable  variance  reduction.  In  general,  we  expect 
that 


Varp.[ZMN  p.(K)]  <<  Varpiid\ZMN  piid(K)]  , 


where 


V^p.[Zmnp.(K)]  =  £ 


_  Aw/(h) 

_  o 

.[s(o/«b  RwvCh) . 

_  Zmn 

which  has  been  verified  to  be  the  case  in  most  of  our  simulations. 
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Probability  Pmn (h),  although  simple  to  compute,  may  result  in  an  unac¬ 
ceptably  high  error  variance.  This  is  typically  the  case  for  GRFs  in  low  tem¬ 
peratures,  as  our  simulations  demonstrate;  therefore,  a  better  choice  for  Pmn  0*) 
is  needed,  which  will  result  in  larger  error  variance  reduction.  This  is  provided 
by  a  theorem,  which  appears  in  [41]. 

Consider  the  following  “entropy”  distance  between  the  Gibbs  distribution 
and  a  probability  distribution  Pmn( h): 


d-iL (kmn ,  Pmn)  -  X 

state  8  h 


KmnO*)  .  |~  rc.w/Ch) 

PWiV(h)  J  ln  [PUN(h) 


PmnQz)  - 


Z  JtAflv(h)ln 

states  h 


^mn  (h) 
Pmn  (h) 


This  is  also  an  Ali-Siluey  type  of  distance,  and  is  minimized  for 
PlWv(h)  =  7tMW(h),  for  all  states  h.  For  the  reasons  stated  above  we  would  like  to 
constrain  the  minimization  problem  and  obtain  a  probability  distribution 
P.vw(h)  which  satisfies 


Pmn (h)  =arg{  min  d,L  (nMN ,  PMN )  )  •  (21) 

Pmn  e  D 

If  we  assume  that  the  Gibbs  distribution  rc^Ch)  has  homogeneous3  LTF’s,  i.e., 


ai;Cc,y,z,co)  =  aix,y,z,  to), 

for  every  (i ,  j )  e  AMN  and  Or,y,z,co)  e  Eh,  then  we  obtain  the  following  theorem. 


Theorem  2:  The  probability  measure  PmnAi)  g  D  with  LTF 


x"(x,y,z,a>) 


31  n  Zmn 

a{x,y,z,  co) 

_  daix  ,y  ,z  ,m)  _ 

z 

ueEH 

31n ZMn 

cyiu  ,y  ,z ,(n) 

da(u,y,z,  to)  ^ 

(22) 


for  every  (x,y,z,to)  e  Eh  and  (i,j)  g  Amn  satisfies  eq.  (21). 

3  This  is  a  reasonable  assumption  in  image  processing  applications  and  simplifies  our 
derivations  considerably. 
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Proof :  The  proof  can  be  found  in  [41]. 


□ 


The  MC-GRF  with  probability  measure  P^N  is  optimal,  achieving  a 
minimum  entropy  distance  from  the  Gibbs  distribution  nMN.  However,  in  gen¬ 
eral,  it  is  suboptimal  with  respect  to  minimizing  the  error  variance,  or 
equivalently  the  distance  dIS(nMN,  Pun )• 

Given  the  LTF  aix  ,y  ,z,w)  of  a  GRF,  our  second  choice  for  the  joint  proba¬ 
bility  distribution  PMnO i)  will  be  a  joint  probability  distribution  PmNQi)  given  by 
eq.  (7),  with  LTF  given  by  eq.  (22).  In  this  case  (see  eq.  (8a)) 


qiJP~(x,y,z,m)  = 


Z 

ueEH  l 


d\nZ, 


MN 


da(u  ,y,z,a>) 


aiu,y,z,t o) 


ainZj 


MN 


<ks(x,y,z,  to) 


(23) 


As  expected,  the  MC-GRF  with  LTF  given  by  eq.  (22)  approximates  the  origi¬ 
nal  GRF  better  than  the  MC-GRF  with  LTF  given  by  eq.  (19),  especially  in  low 
temperatures.  This  may  result  in  significant  variance  reduction;  i.e.,  we  expect 


Varp-[ZMN  p..(K)}  «  Varp.[ZMN  p.(K) ]  , 


where 


_  Aun( h) 

—  O 

.  [state*  h  PmnO*) 

-  %mn 

Varp~[ZMN  p~(K)]  ^ 


which  has  been  verified  by  our  simulation  experiments  for  a  variety  of  GRF’s. 

To  summarize,  we  propose  two  methods  for  calculating  the  partition  func¬ 
tion  ZMn  '■ 


Method  1:  ZMn  =  Ep. 


Amn( h) 


PmnW 


ZMN,P'  - 


V  Z  QmN.P 'fl1*) 
A  *  =  1 


where  h*  ~  P'mn  (h*) 
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Method  2: 


Awiv(h) 

'  1  K  „ 

.PmvOO. 

=>  z  -  =  lim 

MN-P  K->+~ 

S  Q mn.p 

A*  =  i 

where  h*  ~  PuN(hk) , 

with  PmNQ i)  and  Pj^N(h)  being  the  joint  probability  distributions  of  two  MC- 
GRF’s  with  LTF’s  given  by  eq.  (19)  and  eq.  (22)  respectively,  and  Qm  p.(h), 
Q.un  />“(W  by  eqs.  (8b),  (20)  and  (8b),  (23)  respectively. 

IV.  COMPUTATIONAL  ALGORITHMS 

After  deciding  for  the  MC-GRF  joint  probability  distribution  to  be  used  in 
the  Monte-Carlo  partition  function  calculation,  the  next  step  is  to  determine 
how  to  implement  this  calculation  in  a  computationally  efficient  way.  In  the 
following  we  shall  denote  both  Pj^vCh)  and  P^wdi)  as  P^N(h),  i.e,  the  superscript 
Y  will  stand  for  *  or  **.  A  simple  approach  is  to  use  the  following  algorithm: 

ALGORITHM  I: 

1.  Draw  K  statistically  independent  realizations  h* ,  k  =  1,  2,  .  .  .  ,  K,  from 
probability  P^Oi),  given  by  eq.  (7),  with  ty  (x  ,y,z,  o>)  =  x?j(x ,  y ,  z ,  co),  given 
by  eq.  (19)  or  (22).  This  can  be  done  lexicographically,  as  it  has  been  dis¬ 
cussed  in  Section  II. 

2.  Compute  QMN  pT(h*),  for  every  realization  h*,  k  -  1,  2, .  . . ,  K,  by  using  eq. 
(8b)  and  one  of  (20),  (23). 

3.  Compute  ZMN  py(K)  by  using  eq.  (10b). 

The  statistical  properties  of  the  estimator  ZMN  py(K)  can  be  easily  derived. 
Indeed,  ZMN  pl{K)  is  a  sum  of  the  i.i.d.  random  variables  QMN  pT(bi),  QMN  pY(h2) 

,  .  .  .  ,  Qun  pY(h*-)  divided  by  K.  Each  of  these  random  variables  has  finite  mean 
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Kwv.pr =  EpytQMN,py(*11' ) ]  =  zmn  and  finite  variance  pr  =  Var^lQ^  pY(h*)] 
=  £  Qmn  pr^1)  P&nO*)  -  Zun  >  0.  From  the  central  limit  theorem  we  have  that 

staUt  h 

ZMN  pl(K)  is  an  asymptotically  normal  random  variable  with 


2mn  py(K)  -  ZMN 
aMN,py 


-»  iV(0, 1)  , 


as  K  -4  -h».  From  the  strong  law  of  large  numbers  we  also  have  that 

Pr  [  lim  ZMN  py(K)  =  ZMN 1  =  1  ; 

i.e.,  our  estimator  converges  to  the  partition  function  with  probability  one. 
Therefore,  it  converges  in  probability;  i.e., 

Pr  [ 1  zmn.pW  I  >  e]  ->  0  ,  (24) 

as  K  ->  +° o ,  for  every  e  >  0.  Hence,  1UN  pJ(K)  is  an  unbiased  and  consistent 
estimator  of  ZMN.  To  get  a  practical  idea  on  the  accuracy  of  such  an  estimator, 
it  is  worthwhile  calculating  the  sample  variance,  defined  by 

^MN.py^^  =  t  ^MN.py^1^  ~  Q\iN,py(E^  ^  •  (25a) 

A  *=i 

where 

§MN,py^}  =  IF  £  Qmv.pt^1*)  •  (25b) 

A  *=i 

This  is  the  usual  estimator  for  the  variance  of  ZMN  pr(K),  and  can  be  effectively 
used  to  decide  how  many  iterations  K  are  needed  in  order  to  achieve  a 
prespecified  accuracy  in  estimating  Zmn- 

The  computational  complexity  of  Algorithm  1  is  0(KMN),  since  it  requires 
the  generation  of  KMN  samples  drawn  from  the  discrete  conditional  probability 
distribution  Pr[  h i7-  I  A,_ hj,  hi_ij_u  hij_x  ];  therefore,  for  a  large  lattice,  we  may 
not  be  able  to  obtain  accurate  estimates  of  ZMN  in  a  reasonable  time.  This 
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problem  may  be  partially  solved  by  considering  realizations  hlf  ,  which 

differ  only  slightly  from  each  other.  In  this  case,  h*+1  and  QMN  pr(h*+1)  may  be 
easily  calculated  from  h*  and  QMN  py(hk)  by  a  simple  update.  This  can  be  done 
by  generating  a  Markov  chain  on  the  RMN  possible  realizations  of  [if]  with 
prespecified  transition  probabilities  Plf-k+1),  given  by 

pX*,*+n  =  Pr  [H|k+1  =  hj  IH*  =  hj, 

from  state  h,,  at  step  k ,  to  state  hp  at  step  k  +  1.  These  transition  probabili¬ 
ties  form  a  sequence  of  RMNxRMN  transition  probability  matrices  PU), 
k  =  1,2,  ... ,  which  should  satisfy  the  following  relationship: 

(  p'1)1  p(A)=  ( piy  ,  (26) 

where  Pr  is  the  RMN-  dimensional  probability  vector  with  elements  P,Y  = 
P^nO^i)  and  t  denotes  transposition.  If  Pr[  Hi  =  h  ]  =  P^N( h),  and  if  eq.  (26)  is 
satisfied  for  every  k  =  1,  2, ... ,  then  every  state  of  our  Markov  chain  will  be  a 
sample  drawn  from  PJ,N( h).  Therefore,  we  have  to  construct  P(k)  such  that  eq. 
(26)  is  satisfied.  However,  we  have  to  keep  in  mind  that,  in  this  case,  the  reali¬ 
zations  ht,  I12,  . .  . ,  are  no  longer  statistically  independent,  and  a  new  analysis 
is  required  for  the  study  of  the  statistical  properties  of  ZMN  pl(K). 

A  simple  way  to  accomplish  the  previous  ideas  is  to  consider  a  sequence 
hi,  hfc  .  . . ,  of  realizations  such  that  h*+1  differs  from  h*  only  at  one  site  sp, 
this  site  being  chosen  randomly  (among  all  possible  MN  sites  in  AMN),  or,  sys¬ 
tematically  (e.g.,  lexicographically).  Furthermore,  the  random  variable  Hk+l(sp ), 
at  site  sp ,  takes  a  value  h  (sp )  in  EH  drawn  from  a  given  probability  distribu¬ 
tion.  When  this  probability  distribution  is  given  by 

Pr  [fl*+i (sp)  =  h(sp)\Hk(s)  =  h(s),s  e  N,p  ]  >  0  , 

the  resulting  algorithm  is  known  as  the  Gibbs  sampler  [19].  Other  algorithms 
of  a  similar  nature  are  described  in  [41]  and  are  known  as  the  Metropolis’, 
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Barker’s  and  Hastings’  algorithms. 

We  have  compared  the  performances  of  all  these  algorithms,  using  ran¬ 
dom,  or  lexicographic,  site  updating.  As  expected  from  the  analysis  in  [41],  the 
Gibbs  sampler  with  lexicographic  site  updating  gave  the  best  results:  it  con¬ 
verged,  on  the  average,  faster  than  the  other  algorithms,  requiring  less  compu¬ 
tational  time,  and  resulted  in  lower  rejection  rates  (higher  percentage  of  suc¬ 
cessful  updates).  This  result  comes  to  no  surprise,  since  the  Gibbs  sampler  is 
an  importance  sampling  based  procedure  [39],  [42],  [43]. 

The  Gibbs  sampler  is  the  sampling  scheme  used  exclusively  in  our  simula¬ 
tions  and  results  in  the  following  algorithm  for  the  estimation  of  ZMN. 


ALGORITHM  II: 

1.  Generate  a  realization  ht  of  the  MC-GRF  with  probability  PHN(.b),  given  by 
eq.  (7),  with  x  ij(x,y,z,  to )  =  xfiix ,  y ,  z ,  co),  given  by  eq.  (19)  or  (22).  This  is 
done  lexicographically,  as  discussed  in  Section  II.  Calculate  QMN  pY(h1), 
given  by  eqs.  (8b)  and  (20)  or  (8b)  and  (23). 

2.  Set  SUMy  =  QFUNl  =  QMN  p,{h0  and  k  =  1. 

3.  Set  h  =  hk(sp),  with  p  =  (k  -  1)  modulo  MN  +  1,  where  su  s*  ...  ,  sMN,  is  a 
lexicographic  (row  by  row)  ordering  of  the  sites  in  AMN.  Denote  site  sp  by 
d,jl 

4.  Draw  a  value  <j>  e  EH  from  probability  dn  /  ( £f=1  dt ) ,  n  =  1,  2, .  . .  ,  R ,  where 


dt  =  xJOfc, h&\ jM-l j-x , hiikj- 1 )  x  x l , h£l j+i , h£}, Jt*i) 

X  ihx.  j+lihi+i,  j+x .  hfkjKi ,  <t»i ,  hi(kl  j )  , 


for  /  =1,2,...,  R . 

5.  If  <j>  =  h  set  QFUNk  +  i  <-  QFUNk  and  go  to  step  6,  otherwise  compute  the 
ratio 


RATIO  -  Qmh  ptQnk+{)l QMN  py(h-k) , 
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by 


RATIO  -  X  Qi+i. 

"  qij^ih&Ki.hPlj+i.hMj'h)  x  qi+l  j  py(hZl J,h, Aft-i , h$l ^ . 

<li+1,J+upA&\j+i  Mk]+1A,h&\j)  x  q..p^ ) 

X  n  .  ..(hJk)  .hW.,.h.hM  A  x  a  .Jh.hM  :.h&  .■  ,  .  hJk)  .) 


where  qtJ  py(x,y  ,z,a)  is  given  by  eq.  (20)  or  (23).  Set  QFUNk+l  <- 
QFUNk  x  RATIO . 

SUMk+1 

6.  Set  SUMk  +  y  «-  SUMk  +  QFUNk+l.  If  k  +  1  =  K,  set  ZMN  py{K)  =  —  -+*  - 

and  stop;  otherwise,  set  hif+1)  =  (j>,  h£?l)  =  h^,  for  every  im,n)  * 
k  <-  k  +  1,  and  go  to  step  3. 


We  would  like  now  to  emphasize  that,  in  the  case  of  Algorithm  II,  con¬ 
secutive  realizations  of  the  Markov  chain  differ  in  at  most  one  site  and,  there¬ 
fore,  they  cannot  “cover”  rapidly  the  entire  state  space  of  [i7  ],  especially  in  the 
case  of  large  lattices.  A  remedy  for  this  problem  would  be  to  consider  L  <  <  MN 
successive  realizations  of  our  algorithm  as  one  state  transition  of  the  Markov 
chain.  In  practice  however,  we  have  noticed  that  this  modification  does  not 
improve  the  convergence  rate  of  Algorithm  II. 

Once  again,  we  are  interested  in  the  statistical  properties  of  the  estimator 
ZMN  pi(K),  resulting  from  Algorithm  II.  In  order  to  simplify  our  presentation, 
we  shall  consider  Algorithm  II  with  random  (instead  of  lexicographic)  site 
updating,  in  which  case,  the  underlying  Markov  chain  is  characterized  by  sta¬ 
tionary  transition  probabilities.  A  more  general  treatment  (which  is  necessary 
when  lexicographic  site  updating  is  considered)  can  be  found  in  [44],  Both 
cases  result  in  the  same  properties  for  the  estimator  ZMN  py(K).  As  in  the  case 
of  Algorithm  I,  the  estimator  ZMN  py(K)  is  the  sum  of  identically  distributed 
random  variables  QMN  ^(h*),  which  are  characterized  by  a  finite  mean 
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p mn  pi  =  Zmn  and  finite,  nonzero  variance  o2N  py.  Obviously,  tMN  pr(K)  is  an 
unbiased  estimator  of  %mn  ■  The  computation  of  the  variance  of  ZMN  pl(K) 


requires  more  effort.  Observe  that 


Var-py  ZMN  py(K )  -Epy  (  „.  Qmn  pt^*)  -  „2  -®pt  (  X (  Q"an  pt^a)  Zmn)\ 

L  •  J  A  *=i  *=i 


X  £'pThQWiV,PT^,*^-'ZAf^^  +  2  Z  Z  •BpThQWAT,pT(bm)-2Mi>/)(QMAripTfbn)-'Z'Afw)]  . 
**•  Jk=l  *  m 


or 

KVarpy[2MNp1(K)l 

=  Epy  (QMNpyQ*i)-ZMN)2  +2YJ^J^-\LEpy\.QUNipy0^y)QMNpy(hi)]-ZMN 
=  Yj  X  ^MiV.P^bm)  Qmn  pr(bn)  [^^m)  &mn  ~  TjfcvGO  PjftN(hn  ) 

states  hms tales  hn  ’  ’  L 

+  2  Z^PtiN&rnnP^-P&NiKV  ,  (27) 

i=l  A 

where  8mn  is  the  Kronecker  delta,  and 

P^>  =  Pr[Hi  =  h„  lH1  =  hw  ]. 


Equation  (27)  cam  be  easily  written  in  a  matrix  form  as 

f  K-l 

KVarpy[ZMN  pTCO  =  ( £Y)‘  PY  -  PYA  +  2  PY  £ 

’  L  /=I 

where  all  the  matrices  are  RMNx  RUN  matrices.  Matrix  PY  is  a  diagonal  matrix 
with  elements  the  RMN  probabilities  PjfrN(h),  h  e  Eft*.  Matrix  A  has  RMN  identi¬ 
cal  rows  which  are  equal  to  the  diagonal  of  PT,  whereas,  P'  is  the  l  power  of 
the  transition  probability  matrix  P.  Finally,  Q1  is  an  RMN-  dimensional  vector, 
with  elements  QMN  pT(h),  h  e  E$N.  For  an  ergodic  transition  probabihty  matrix 
P  the  fundamental  matrix  F,  given  by  F  =  (I  -  P  +  A)"1,  exists,  and  [45] 


^(P'-A)  Qf, 
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F  =  I  +  lim  *£  (P'  -  A) .  (29) 

K  — *  +°*  [s\  ^ 

From  eqs.  (28)  and  (29)  we  obtain 

\im_KVarpl\2MN  py(K^  =  (&)*  [2P*F  -  P*-  PTa]  gr  , 
or 

lim  Varpy  \zuN  py(K )  =  0  .  (30) 

From  eq.  (30)  and  the  Tchebycheff  inequality  we  obtain  eq.  (24);  therefore, 
Algorithm  II  results4  in  an  unbiased  and  consistent  estimator  for  ZMN  . 

We  would  also  like  to  get  a  practical  idea  of  the  accuracy  of  the  resulting 
estimate.  However,  since  the  samples  h*  are  no  longer  independent,  we  cannot 
use  the  sample  variance  given  by  eq.  (25).  In  the  case  of  Algorithm  II,  eq.  (25) 
will  most  probably  underestimate  the  variance,  since  it  does  not  take  into  con¬ 
sideration  the  correlations  among  consecutive  realizations.  To  overcome  this 
difficulty,  we  may  use  the  ideas  in  [48]. 

In  many  instances  ZMN  ->  -h»  with  increasing  lattice  size  (i.e.,  as 
M,N-*+°°).  Additionally,  for  many  realizations  h,  P^N(h)  0,  whereas, 
Amn0 b)  — >  -k»  as  M ,  N  -»+«..  Therefore,  eq.  (5c),  or  equivalently,  eq.  (11),  may  be 
violated.  Since  most  applications  require  use  of  rather  large  lattices,  we  have 
to  examine  the  problems  introduced  by  this  behavior  and  modify  our  Monte- 
Carlo  calculation  procedure.  It  is  clear  that,  in  these  cases,  a  simple  partition 
function  calculation  algorithm  (like  Algorithm  I,  or  Algorithm  II)  will  suffer 
from  overflow  problems.  A  first  step  towards  the  remedy  of  this  problem  is  to 
estimate  the  quantity 

f(ZMN)  =  \n(ZMN) ,  (31) 

4  The  central  limit  theorem,  developed  in  [46],  [47],  is  directly  applicable  in  our  case,  and, 
therefore,  our  estimator  is  also  asymptotically  normal. 
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instead  of  ZMN,  by  (see  also  eqs.  (8b)  and  (10b)) 


/ Ktpy(ZMN)  -  In  (if ))  ln(Q^^  pT(hi))  +  In 


K  Q 


MN.py 


(hk) 


K e^.prOn) 


,  (32a) 


with 


In  (Qmn  p1<hd)  =  S  2  In  [9  pT(^«  A/iJ.  Jf  •  A/Zi-i >1  ■  (32b) 

’  i=l 7=1  1  J 


where  q{j  pr(x,y,z, co)  is  given  by  eq.  (20)  or  (23).  This  is  a  reasonable  approach 
for  two  main  reasons:  (a)  f  (ZMN)  will  not  necessarily  become  infinite  as 
M  ,N ,  ZMN  -»+<*>  (see  for  example  the  case  of  the  Ising  model,  discussed  in  Sec¬ 
tion  VI),  or,  at  least,  it  will  be  much  smaller  than  ZMN ;  and,  (b)  we  are  usually 
interested  in  computing  the  logarithm  of  the  partition  function  rather  than  the 
partition  function  itself  (e.g.,  in  the  case  of  maximum  likelihood  parameter  esti¬ 
mation).  However,  computing  the  summation  in  eq.  (32a)  may  still  result  in 
overflow. 

To  overcome  this  problem  observe  that  (see  also  eqs.  (8b),  (9)  and  (31)) 


f  (%mn  )  - 


I  Q^,Pr(b)P^(h) 

states  h 


MN 


In 


states  h 


M  N 

nn 

»=U=1 


j  ’  ^«~1.  j-U^i  ,  j-0 

qij,Py 


n  (max) 

qij,P'1 


PM h) 


MN 


In 


z  h) 

states  h 


=  In  ^ 


MN 


In 


Z  Qww.pr(h)^(h) 

states  h 


(33) 


where 


n  (max) 

qU,Py 


=  max 
Cx ,*,*,<•)«  E$ 


( x,y,z,a> ) 


qmax  ,Pf 


M  N 

nn«S' 

i.ij.i 


i 

MN 


page  24 


and 


QmN.PT^  -  nn  (mat) 

■=ly=l  Qij'P  T 


From  eq.  (33)  we  see  that,  in  order  to  compute  f  (ZUN)  we  have  to  calculate  the 
expectation  Epy[QMN  M(h)].  This  can  be  done  by  a  simple  modification  of  the 
previous  two  algorithms.  Notice  that,  in  this  case,  Monte-Carlo  calculations  can 
proceed  with  no  overflow,  even  in  the  case  of  large  lattices,  because  Q^N  pT(h) 
<  1<-h>°,  for  all  realizations  h;  therefore,  condition  (11)  is  satisfied.  Equation 
(33)  yields  the  following  Monte-Carlo  estimate  for  f{ZMN) 


f K.PT^Mn)  -  fyipf^ZMN.PT^))  -  ^Qnax.pT  + 


MN 


In 


K  £  Qmn.pt 

A  k  =  i 


(h*) 


where  h*  are  samples  drawn  from  the  joint  probability  distribution  W h). 

It  is  worthwhile  noticing  that  the  estimator  In  (ZMN  py(K))  will  only  asymp¬ 
totically  (i.e.,  for  K  -»+«)  yield  an  unbiased  estimate  for  In (ZMN)  [49].To  see 
this,  and  for  the  case  of  Algorithm  I,  use  the  central  limit  theorem  for 
ZMN  py(K)  and  the  fact  that  the  function  ln(- )  has  a  nonzero  derivative  at  ZMN 
to  prove  that  \n(ZMN  py(K))  is  asymptotically  normal  with  mean  In (ZMN)  and 
variance  UK  (aMN  pylZUN)2,  where 

<*MN.pys  X  Q*N'Py(h)P,«N(h)-Z$N  . 

stales  h 


A  similar  result  can  be  proved  for  the  case  of  Algorithm  II.  Finally,  we  can 
easily  show  that  UMN  \n{2,^N  pfK))  is  a  consistent  estimator  for  I /MN  m(Z.WjV  ), 
when  either  one  of  the  two  algorithms  is  used. 

The  practical  implementation  of  the  previous  ideas  requires  the  develop¬ 
ment  of  fast  algorithms.  There  are  many  ways  to  improve  performance.  For 
example,  much  time  is  devoted  to  generating  uniformly  distributed  random 
numbers.  The  development  of  a  good,  fast  and,  probably,  machine-dependent 
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random  number  generator  will  definitely  result  in  computational  savings. 
Another  problem  is  the  fact  that  the  rejection  rate  approaches  100%  in  low 
temperatures.  To  avoid  “being  stuck”  at  the  same  realization  for  a  long  time, 
one  might  try  to  use  some  “a-priori”  knowledge  for  the  MC-GRF  parameters 
and  perform  various  simulations  with  different  suitable  configurations  of  the 
state  space,  averaging  the  obtained  results  with  appropriate  weights.  An 
ingenious  idea  to  overcome  a  similar  problem  is  the  “n-fold  algorithm”  proposed 
in  [50].  Since  in  low  temperatures  most  GRF’s  will  favor  realizations  which  are 
characterized  by  large  clusters,  one  can  try  to  break  these  clusters  into  smaller 
ones  and  then  assign  a  random  spin  value  to  the  whole  cluster.  This  idea  is 
used  in  [51]  to  obtain  an  efficient  and  fast  Monte-Carlo  simulation  algorithm. 
It  is  also  important  to  notice  that  Algorithm  II  can  be  directly  implemented  on 
the  Monte-Carlo  simulation  machines  developed  in  [52-55].  These  are  spf  rial 
purpose  computers  developed  for  the  study  of  various  properties  of  GRF’s. 
Implementation  of  our  algorithm  on  such  machines  ma;/  result  in  a  substantial 
reduction  of  computational  time. 

Traditionally,  stochastic  methods  are  used  to  simulate  GRF’s  and  calcu¬ 
late  various  quantities  related  to  them.  Our  method  fits  to  this  framework. 
However,  some  researchers  have  adopted  deterministic,  or  pseudo- 
deterministic,  evolution  algorithms,  like  the  microcanonical  simulation  algo¬ 
rithm  in  [56],  [57],  For  such  algorithms,  it  is  much  more  difficult  to  prove  con¬ 
vergence  and  study  statistical  properties  of  the  resulting  estimators.  Addition¬ 
ally,  it  is  not  clear  if  they  can  result  in  significant  computational  savings.  How¬ 
ever,  we  believe  that  it  is  extremely  interesting  to  pursue  many  of  these 
methods  further.  We  haven’t  done  so  here  though,  since  our  major  objective  is 
to  present  fundamental  ideas,  and  a  simple  algorithm,  for  the  calculation  of 
Amn  ■ 
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V.  CALCULATION  OF  PARTITION  FUNCTION  DERIVATIVES 


So  far,  we  have  discussed  the  problem  of  computing  the  partition  function 
ZMN  of  a  GRF.  We  have  proposed  two  methods  which  require  sampling  from 
two  different  MC-GRFs,  namely,  P^s  and  Pmn-  Method  2  requires  the  computa¬ 
tion  of  the  R 4  first-order  derivatives  of  the  logarithm  of  the  partition  function 
with  respect  to  the  LTF  (see  eq.  (22));  i.e.,  it  requires  computation  of  the  vector 


=  ( 


3a(x,y,z,o>)  ’ 


(x,y,z,a)  e  Eh  }  . 


(34) 


The  analytical  computation  of  eq.  (34)  is  not  possible  in  general  (except  in  some 
special  cases).  However,  this  calculation  can  be  approximated  via  Monte-Carlo 
simulations  by  using  a  scheme  similar  to  that  used  in  statistical  mechanics  for 
the  calculation  of  various  thermodynamic  properties  of  large-scale  systems  [58]. 
In  addition  to  the  derivatives  in  eq.  (34)  we  may  also  want  to  compute  the 
internal  energy  (which  is  related  to  first-order  derivatives  of  the  partition  func¬ 
tion)  and  the  specific  heat  (which  is  related  to  second-order  derivatives  of  the 
partition  function)  of  a  general  GRF.  These  quantities  provide  valuable  means 
for  detecting  and  studying  phase  transitions.  All  these  computations  can  be 
approximated  by  a  simple  Monte-Carlo  simulation  algorithm,  as  it  will  be 
demonstrated  in  this  section. 

Let  us  study  the  calculation  of  the  derivatives  of  the  logarithm5  of  the 
partition  function  ZMN  with  respect  to  the  LTF  aix,y,z, co).  We  shall  limit  our 
analysis  to  the  computation  of  first-  and  second-order  derivatives.  However, 
our  presentation  can  be  trivially  extended  to  include  the  computation  of  any 
higher-order  derivative.  In  the  following,  we  shall  denote  by  v*(x,y,2,<o)  the 

5  We  prefer  to  calculate  the  derivatives  of  the  logarithm  of  the  partition  function,  instead 
of  the  derivatives  of  the  partition  function  itself,  for  reasons  similar  to  those  explained  at  the 
end  of  Section  IV. 
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total  number  of  elementary  squares  (hif\  hf^  j,  h^j.i ,  )  =  (x,y,z, to) 

which  appear  in  a  realization  h*  of  the  GRF  [H],  From  eqs.  (2b)  and  (3a)  we 
have  that 


ZMN  =  I  FI  oiq,u,v,w)v‘(q'u’u’w)=  ZAMNQii), 

statesb'  (q iU  ,u ,w)<=  Efi 


(35a) 


states  h 


where 


AMN(hi)=  n  c{q,u,v  ,w)''l{q'u'u,w) . 

(q  ,u  ,u  ,w)e  Efh 

Differentiating  eq.  (35a)  with  respect  to  the  LTF  dx  ,y  ,z  ,co)  we  obtain 


Z^ix,y  ,z,a>) 


1  dlnZi fN 

MN  da(x  ,y  ,z  ,co) 


Z  vlOc,y,z,ai)AMN(lil) 

^  states  h{ 

MN  dx ,y ,z ,m)  Z  Awiv(hi) 

states  hj 


(35b) 


_1 _ 1 

MN  dx,y,z, to) 


•  z 

states 


vt(x,y,z,  co)  3tMN(hz)  , 


(36) 


whereas,  differentiating  eq.  (36)  with  respect  to  dq  ,  u ,  v ,  w ),  we  obtain  (see  also 
eq.  (35)) 


Ztfh(x,y,z,<a;q,u,v,w) 


1 _ d2ln Zmn _ 

MN  d  dx  ,y  ,z  ,co)  d  dq  ,u,u  ,w) 


_1 _ 1 _ 

MN  dx  ,y  ,z  ,(o)dq  ,u  ,v  ,w) 


x 


X 


VzU.y.z.co)  vt(q,u,v,w)  -  \>i(x  ,y  ,z  ,&>) 


Kmn  (l1: ) 


Z  ^i(q,u,v,w)nUN(hi) 


states  h, 


Z  v/(*,y,z,e>)7rWN(hz) 

states  h( 


L,  (37) 
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where 


jl  ,  \{{x,y,2,a>)  =  (.q,u,v,w) 
[0  ,  if  (x,y,z,co)  *  (q,u,v  ,w) ' 


The  quantities  defined  by  equations  (36)  and  (37)  can  be  efficiently  calcu¬ 
lated  via  Monte-Carlo  simulations  based  on  drawing  samples  directly  from  the 
Gibbs  distribution  (2).  This  is  a  quite  difficult  problem  in  practice,  whose  solu¬ 
tion  can  only  be  obtained  asymptotically,  therefore  the  resulting  Monte-Carlo 
estimates  are  biased.  In  practice,  we  will  generate  a  Markov  Chain  of  realiza¬ 
tions  h* ,  k  =  1,2 with  suitable  transition  probabilities,  that  will 
asymptotically  be  distributed  according  to  the  Gibbs  distribution,  i.e., 
h*  ~  *mn  ,  or  lim  Pr[  H  =  h*  ]  =  nUN(hk),  given  by  eq.  (2).  In  order  to  reduce 

*  ~  *  -*4— 

the  variance  of  the  resulting  estimates,  we  discard  the  Kx  initial  configurations 
of  this  Markov  Chain,  where  Kx  is  a  sufficiently  large  number,  such  that 
Pr[  H  =  h*  ]  =  KMN(hk),  for  k  >  Kx.  Then,  we  consider  the  K2  next  states  of  the 

Markov  Chain,  i.e.,  h*,  k  =  Kx+l,  Kx  +  2 . Kx+K2,  on  which  we  form  ergodic 

averages.  This  results  in 


%mn(x  ,y  ,2,“;  Ki,k2)  - 


_i _ l  l 

MN  aix,y,z,a>)  K2 


Kx+K2 

Z  v*(x,y,z,w)  , 

*  =/Cj  +  l 


(38a) 


and 


Z$n(x  ,y  ,z  ,co;  q  ,u  ,v  ,w ;  KUK2)  = 


MN  a{x ,y  ,z  ,(o)dq  ,u  ,v  ,w) 


^  Kl  +  K2  r 

x  \  —  X  v*  (x  ,y  ,z ,  <o)  v* (q ,  u ,  w ,  u; )  -  v*  (x  ,y , z ,  to)  83^"" 


2  *=x1  +  iL 


_  _ 

r  i 

•\ 

Kt+K2 

Ki+K2 

Z  v*(x,y,2,  to) 

Z  vk{q,u,v,w) 

+ 1 

k=Kx*  1 

■  J 

.  J 

J 

(38b) 
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which  are  the  Monte-Carlo  estimates  of  the  the  first-  and  second-order  deriva¬ 
tives  of  the  logarithm  of  the  partition  function  with  respect  to  the  LTF’s. 

There  are  many  important  applications  which  require  the  calculation  of 
the  partition  function  and  its  derivatives.  An  anticipated  application  is  the 
development  of  an  optimal  algorithm  for  the  maximum-likelihood  estimation  of 
the  LTF  of  a  GRF  from  a  given  realization.  We  shall  not  expand  on  this  subject 
here;  instead,  we  shall  illustrate  the  use  of  these  calculations  for  the  computa¬ 
tional  study  of  phase  transitions. 

The  internal  energy  EMN(T)  and  the  specific  heat  CMN(T),  given  by 


Emn(T ) = 


1  „2  dlnZarv 
MN  ZT 


(39a) 


and 


<WD  = 


3Eun(T) 

ZT 


-  -^Emn(T)  + 


1  2  d2ln ZMN 

mn  ar2 


(39b) 


respectively,  are  two  important  thermodynamic  quantities,  associated  with  a 
GRF,  which  allow  the  study  of  phase  transitions  (5].  If  Tc  is  a  temperature 
such  that 


lim  lim  CMN(T)  =  +°°,  (40) 

T  -*TC  M,  N  -*•*- 


then  we  say  that  the  GRF  is  in  phase  transition  at  critical  temperature  Tc .  The 
study  of  phase  transition  is  very  important,  being  the  subject  of  many  discip¬ 
lines,  including  statistical  mechanics  [2],  [3-5],  [29],  [34],  [44],  [59],  probability 
and  information  theory  [2],  [32-34],  [44],  [60],  [61]  and  image  processing  [12], 
[34],  [44].  From  eqs.  (39)  and  (40)  we  see  that  the  critical  temperature  Tc  is  the 
temperature  at  which  the  internal  energy  EMN(T),  or  the  first-order  derivative 
d\nZMN/dT ,  is  discontinuous.  Observe  that 


ovu^MN 


( x,y,z,a)e  Efi 


idc(x,y  ,z,co)  3T 
1 1 


dT 


(41) 


page  30 


and,  therefore,  the  internal  energy  EMN(T)  can  be  approximated  by  (see  also 
eqs.  (36),  (38a),  (39a)  and  (41)) 


£UN(T-,KuK2)  =  T2  £  tj&(x,y,z,*\KuKd 

(x  .y  .  o»'e  Ef} 


da(x,y,z,  co) 
dT 


_  rp2 


_L _ 1_ 

MN  K2 


K\+K2 


I 

k=Kl  + 1 


X  v*(x,y,z,co) 

(x,ji,z  ,u)e  Efi 


d\na(x  ,y  ,z  ,(o) 
dT 


(42) 


At  the  critical  temperature,  the  second-order  derivative  d2\nZMN/dT2  will  neces¬ 
sarily  become  infinite.  From  eqs.  (36)-(38),  (39b)  and  (42)  we  can  easily  show 
that  the  specific  heat  CMN(T)  can  be  approximated  by 

Cmn{T-,kx,k2)  =  t2  £  l2tikx,y,Z,co;KuK2) 

( x,y  ,z  ,<o)  e  Efj  l 


2  da(x,y,z,co)  d2<y(x, y,z,  co) 
T  dT  dT2 


do(x,y,z,(o) 

dT 


X  2mn0c ,y  ,z,<t>;q ,u,v  ,w,  KuK2)  3- ’ U  ’ W- 

(q ,u,  u  ,w  )e  Efi 


1  1  1 
T2  MN  K2  1 


x1+x2 

X 


(x,y,z,a)e  E$ 


1 

^2 


^  ^,2  31no(x,y,z,co) 
X  X  v*(x,y,z,co)  7^ - 1 — 

+  l  (.x,y,z,u)z  Efl 


(43) 


To  effectively  compute  the  Monte-Carlo  estimates  given  by  eqs.  (38a),  (42) 
and  (43)  we  employ  the  Gibbs  Sampler  algorithm  with  a  lexicographic  site 
updating,  by  generating  a  Markov  Chain  which  asymptotically  reaches  the 
Gibbs  distribution.  The  overall  procedure  is  summarized  in  the  following  Algo¬ 
rithm  III. 
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ALGORITHM  III: 


Initialization 

1.  Generate,  lexicographically,  a  realization  ht  of  the  MC-GRF  with  probabil¬ 
ity  PmnQA),  given  by  eq.  (7),  with  xu(x ,y,z,  co)  =  x*(x ,y,z,  co),  given  by  eq. 
(19).  Set  k  <-  1. 

2.  Set  A  =  hk(sp),  wdth  p  =  (k  -  1)  modulo  MN  +  1,  where  S],  .  .  .  ,  is  a 

lexicographic  ordering  of  the  sites  in  AMN.  Denote  site  sp  by  (i ,  j ). 

3.  Draw  a  value  <$>  from  probability  dn  /  ( £f=1  di) ,  n  =  1,2, ...  ,  R,  where 


dt  =  aitfi ,  h^l  j ,  h^l  j-i ,  )  x 

x  oihtfl  , k£l  j-! )  x  y+i ,  V.% » 4>i .  j ) . 


for  I  =1,2,...,  i? . 

4.  Set  htf+l)  =  h£‘n+1)  =  h™,  for  every  (m ,  n)  *  (i ,  j).  Set  k  <-  k  +  1. 

5.  If  A  =  +  2  compute  vXl+1Cr,y  ,2,00),  for  all  (x,y  ,z,co)  e  and 


TERMKi  +  l  = 


T2  X  vx1+i(*,y,z,“) 


dInq(;ic,y,z,co) 

dT 


(44) 


Then, 

derKx+l(x,y,z,m)  =  vKl+l(x,y,z,  co),  for  every  (x  ,y  ,z  ,w)  e  E$  , 
SUM  1*I  +  1  =  TERMKl+l  ,  SUM2Kl+1  =  TERM^  +  1  ; 
otherwise,  go  to  step  2. 


Main  Iteration 

6.  Set  A  =  A*(sp),  with  p  =  (is  -  1)  modulo  MAT  +  1,  where  si,  s*  ...  ,  sMN,  is  a 
lexicographic  ordering  of  the  sites  in  AMN .  Denote  site  sp  by  ( i ,  j ). 

7.  Draw  a  value  <t»  from  probability  dn  /(£f=i  dt) ,  n  =  1,  2 , ....  R,  where  dt, 
l  =  1,  2, .  . . ,  R  are  given  by  eq.  (44). 
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8.  If  0  =  hif } ,  go  to  step  9;  otherwise,  compute  vk  +  1(x,y  ,z,co)  of  the  new  reali¬ 
zation,  for  all  {x ,y ,z ,<a)  e  Eh-  Also,  compute: 

TERMk  +  1  =  T2  £ 

(x,y  ,z, a)  e  Efi 

9.  Set 

derk  +  l(x  ,y,z,(a)  <-  derk(x  ,y  ,z ,<o)  +  v*  +  1(x  ,y  ,z ,  a) ,  for  every  (x,y,z,a)  e  , 
St/M  1*  +i  <-  SUM  1*  +  TERM*  +1  ,  SUM2k  +  l  <-  SUM2k  +  TERMk+  k  . 

10.  If*  +  1=K1  +  K2,  set 


2&t,(x,y,z,*;KuKd  = 


MN  aix,y,z,o. >)  K2 


1  -]£-derKi+K2(x,y,z,a), 


Em(T;  *,.*„)  =  £  SMI,,.,, , 


Cmn(T\  K\,K2)  - 


r2  mn  if2 


SUM2/c1+k2  ~  SUM1k1+k2 


and  stop;  otherwise,  set  hif+l)  =  0,  *^*+1)  =  h**n\  for  every  {m,n)  *(i,  j),  k 
k  +  1,  and  go  to  step  6. 


The  obtained  estimates  will  be  consistent,  but  only  asymptotically  unbiased. 


VI.  SIMULATION  EXPERIMENTS 

We  have  demonstrated  the  fact  that  the  partition  function  of  a  general 
GRF  can  be  calculated  by  using  the  Monte-Carlo  schemes  discussed  in  Section 
IV.  Two  major  questions  have  to  be  answered  at  this  point.  The  first  is  how 
well  the  proposed  algorithms  work.  This  question  is  essential,  because,  in  prac¬ 
tice,  there  are  many  factors  that  can  prevent  us  from  getting  the  expected 
results.  Some  of  these  factors  are:  (a)  the  quality  of  the  random  number  gen¬ 
erator  used  in  the  program,  (b)  the  truncation  errors  incorporated  in  the 
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computation  of  various  sums;  (c)  the  number  K  of  samples  used  for  the  compu¬ 
tation  of  the  estimates;  and,  (d)  the  numbers  Kx  and  K2  used  in  Method  2  for 
the  calculation  of  the  first-order  derivatives  of  the  partition  function  (see  Sec¬ 
tion  V).  The  mishandling  of  some,  or  all,  of  these  factors  may  give  misleading 
results.  The  second  question  concerns  the  relative  merits  of  Method  1  and 
Method  2.  We  are  interested  in  knowing  which  method  gives  more  accurate 
results,  and  which  is  more  computationally  efficient.  Additionally,  we  would 
like  to  demonstrate  the  fact  that  sampling  from  the  probability  distribution  of 
a  MC-GRF  with  the  i.i.d.  choice  for  the  LTF  fails  to  give  reasonable  estimates, 
whereas,  our  methods  constitute  a  vast  improvement  over  the  i.i.d.  case.  All 
these  can  be  accomplished  by  comparing  our  computational  results  with 
analytical  results  obtained  by  either  performing  the  required  summations  over 
all  Rmn  possible  realizations  of  the  GRF  [H  ],  or  by  using  analytically  known 
solutions.  The  first  approach  is  limited  to  GRFs  defined  over  small  lattices, 
whereas,  the  second  approach  is  usually  limited  to  the  2-D  Ising  model  with  no 
external  magnetic  field  and  nearest  neighbor  interactions  [2],  [5],  [32-34],  [59], 
[61-66].  We  shall  present  various  comparisons,  by  employing  both  approaches, 
and  demonstrate  that  the  proposed  methods  provide  highly  accurate  results. 

In  the  first  set  of  experiments  we  consider  five  binary  GRF’s,  at  different 
temperatures,  with  a  second-order  neighborhood  system  and  a  homogeneous 
LTF,  defined  initially  over  a  small  rectangular  lattice  of  4x4  sites.  These 
GRFs  are  fully  described  in  terms  of  16  parameters  a  =  {  a(ac,y,z,to),  x,  y ,  z,  co 
=  0, 1  }.  A  program  has  been  written  which  estimates  the  partition  function  by 
using  Algorithm  II  and  by  sampling  from  the  three  different  probability  distri¬ 
butions  P$n,  P'mn »  and  P'mn ■  Estimates  of  V2ln ZMN,  EMN(T)  and  CWJV(r)  are  also 
obtained  by  using  Algorithm  III.  All  estimated  quantities  are  compared  with 
the  actual  ones,  which  are  calculated  by  performing  summations  over  all  possi¬ 
ble  64,000  states.  Various  simulation  experiments  have  been  performed  for  the 
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three  different  sampling  probability  distributions  and  for  a  wide  range  of  tem¬ 
peratures.  The  results  are  extremely  accurate  and  demonstrate  the  fact  that 
both  our  methods  work  well,  converging  to  the  correct  values  within  a  small 
number  of  iterations  (usually  less  than  a  thousand).  The  method  based  on  the 
i.i.d.  choice  fails,  especially  at  low  temperatures.  The  number  of  iterations 
greatly  depends  on  the  choice  of  the  LTF  of  the  underlying  sampling  probabil¬ 
ity  distribution,  as  well  as  on  the  temperature.  At  high  temperatures,  our 
methods  converge  extremely  fast,  whereas,  at  low  temperatures  many  itera¬ 
tions  are  necessary.  This  is  a  predictable  behavior,  since,  at  high  temperatures, 
the  GRF  and  the  MC-GRF  are  both  “close”  to  an  i.i.d.  random  field,  whereas, 
at  lower  temperatures,  the  MC-GRF  is  only  an  approximation  of  the  original 
GRF  [38].  This  is  especially  obvious  when  the  GRF  exhibits  strong  diagonal 
interactions  between  sites  (i,j)  and  (i~l,j+l),  in  which  case.  Method  1  does  not 
perform  much  better  than  the  i.i.d-  case;  however,  Method  2  is  always  superior. 
This  can  be  seen  by  calculating  the  error  variances  of  the  different  sampling 
schemes  by  performing  the  necessary  summations  over  all  the  64,000  states, 
and  by  assuming  that  Algorithm  I  is  used.  These  calculations  show  that 
Method  2  results  in  the  lowest  variance,  whereas,  in  most  cases,  Method  1  has 
lower  variance  than,  the  i.i.d.  case. 

Figure  1  depicts  realizations  of  the  five  GRFs  considered  here.  Each  of  the 
five  rows  corresponds  to  realizations  of  one  128  x  128  site  GRF  at  five  different 
temperatures.  Phase  transition  is  apparent  in  these  realizations.  As  the  tem¬ 
perature  drops  from  high  ( T  >  Tc ;  realizations  look  random)  to  low  (T  <  Tc ; 
realizations  look  ordered  and  structured),  the  qualitative  behavior  of  the  ran¬ 
dom  field  changes,  and  short  range  interactions  among  sites  develop  into  long 
range  ones.  Table  I  depicts  the  different  temperatures  used  for  the  realizations 
of  Fig.  1,  whereas,  Table  II  depicts  the  values  of  vector  a.  In  Figs.  2A-2E  a 
comparison  of  the  exact  value  of  f(Z 4,4)  and  the  estimated  one,  obtained  by 
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using  Algorithm  II  with  all  three  sampling  methods  (i.i.d.,  Method  1,  and 
Method  2),  is  depicted.  In  this  case  K  =  Kx  =  K2  =  1,600.  In  many  instances,  the 
i.i.d.  case  fails  to  give  a  reasonable  result,  whereas,  our  methods  give  perfect 
results.  The  exact  internal  energy  E^T)  and  specific  heat  Citi(T)  are  also  dep¬ 
icted,  as  well  as  their  corresponding  estimates,  obtained  by  using  Algorithm 
III.  A  high  degree  of  accuracy  is  achieved.  Figures  3A-3E  illustrate  the 
behavior  of  our  methods  as  compared  to  the  i.i.d.  case.  The  minimum  number 
of  iterations  KP  necessary  to  achieve  a  certain  degree  of  confidence  in  the  qual¬ 
ity  of  our  estimates  (a  5%  confidence  interval),  is  plotted,  in  a  logarithmic 
scale,  as  a  function  of  temperature  T,  for  the  five  different  GRF’s  of  Table  II.  If 
we  assume  an  i.i.d.  sampling  scheme  (like  that  of  Algorithm  I)  and  that  an 
asymptotically  normal  distribution  is  achieved,  then  we  can  easily  show  that 

VarP  (Qmv  p  ) 

KP  =  int  100.0 - =>  Pr[ZMN:P(K)  e  (0.8  ZMN,  1.2  ZMN)  ]  >  0.95  , 

ZMN  J 

for  K  >Kp.  Method  2  results  in  substantial  improvement  (with  respect  to  the 
minimum  number  of  iterations)  over  the  i.i.d.  case. 

From  the  first  set  of  experiments,  we  can  conjecture  that  Method  2  is 
superior  to  Method  1.  Now,  we  would  like  to  see  whether  this  is  true  as  the  lat¬ 
tice  size  grows.  In  the  sequel,  we  consider  the  same  five  GRFs  on  larger  lat¬ 
tices  (up  to  128x128  sites).  We  use  Algorithms  II  and  III  to  obtain  estimates  of 
Zmn(K)  (sampling  from  Pjjft,  P'MN  and  P'mn),  of  EMN(T),  and  CMN(T).  Although 
exact  analytical  results  are  unavailable  for  general  GRFs,  we  can  draw  some 
conclusions  about  the  efficiency  of  our  two  methods  by  checking  the  behavior  of 
the  estimated  curves  of  f(2MNiP)  versus  the  temperature  T  as  the  lattice  size 
becomes  larger.  Such  curves  are  depicted  in  Figs.  4A-4E  (for  16x16  sites),  Figs. 
5A-5E  (for  32x32  sites),  Figs.  6A-6E  (for  64x64  sites),  and  Figs.  7A-7E  (for 
128  x  128  sites),  for  the  five  GRF  models  of  Table  II.  In  all  cases 
K  =  Kx  =  K2=  10,000AflV .  The  shape  of  these  curves  indicates  that  Method  1 
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breaks  down  at  low  temperatures  (for  T  <TC),  being  comparable  to  Method  2  at 
high  temperatures  (for  T  >  Tc ).  However,  it  is  impossible  to  judge  the  accuracy 
of  Method  2,  especially  at  low  temperatures  and  at  temperatures  around  the 
critical  temperature,  since  we  have  no  exact  results  to  compare  with.  It  seems 
though,  that  Method  2  gives  reasonably  good  estimates,  except  in  the  case 
when  the  random  field  exhibits  strong  diagonal  interactions  between  sites  (i,j) 
and  (i -1,_/+1).  This  is  evident  from  Figs.  1,  4E,  5E,  6E,  and  7E;  the  estimated 
curves  of  the  partition  function  versus  T,  obtained  from  our  simulations,  are 
not  “smooth  enough”  at  low  temperatures,  implying  large  error  variance.  Notice 
however  that,  the  remarkably  good  results  obtained  from  the  second  set  of  our 
experiments,  described  next,  are  very  encouraging  and  strongly  favor  Method  2 
over  Method  1.  Finally,  to  obtain  a  rough  idea  about  the  location  of  the  critical 
temperature  Tc  of  the  GRF’s  under  consideration,  we  plot  the  internal  energy 
and  specific  heat  curves  of  the  GRF  models,  for  different  lattice  sizes  (8x8, 
16x16,  and  32x32).  The  resulting  graphs  are  depicted  in  Figs.  8A-8E;  they 
allow  us  to  conclude  that  the  break  down  point  of  Method  1  is  indeed  very  close 
to  Tc.  The  critical  temperature  corresponds  to  the  peak  of  the  specific  heat 
curve  as  the  lattice  size  grows  (see  eq.  (40)). 

In  the  second  set  of  our  experiments,  we  compare  our  results  with  analyti¬ 
cally  known  solutions.  We  consider  a  special  case  of  a  GRF,  the  two- 
dimensional  Ising  model  with  no  external  magnetic  field  and  nearest  neighbor 
interactions,  defined  on  a  rectangular  lattice.  This  model  has  been  first  intro¬ 
duced  by  the  German  physicist  Ernst  Ising  in  his  attempt  to  statistically  for¬ 
mulate  the  phenomenon  of  ferromagnetism.  Although  the  Ising  model  has  a 
very  simple  LTF,  it  has  been  enjoying  considerable  attention,  because  it  is  one 
of  the  rare  non-trivial  (i.e.,  it  exhibits  phase  transition)  GRF  models  whose  par¬ 
tition  function  and  its  derivatives  can  be  computed  analytically.  The  neighbor¬ 
hood  system  for  such  a  model  is  the  first-order  system  N^,  defined  in  Section 
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II,  whereas,  R  =  2,  EH  =  (-1,  +  1}  and6 


a(X,  Y ,Z,£l)  =  exp 


y(BXY  +AXQ) 


(45) 


If  we  assume  torodial  boundary  conditions,  we  can  compute  the  partition  func¬ 
tion  analytically,  as  it  has  been  done  initially  by  Onsager  [59].  If  we  define 
parameters  a=AlT  and  $  =  BlT  ,  then  the  partition  function  of  the  Ising  model 
is  given  by  [32],  [63]  (for  all  temperatures  T  such  that  T  *TC) 

MN_  r  N  N 

Zmn  =  -5-  (2  sinh2oc)  2  J  [  2  cosh(  —  y2r  )  ]  +  [  2  sinh(  — y2r  )  1 

J  [r=l  i  r= 1  Z 

+  PJ  [  2  cosh(  — - y2r_i )  ]  +  n  [  2  sinh(  — y2r_1 )  1  k  (46a) 

r=l  A  r-l  “ 


where 

coshy,-  =  cosh2a*  cosh2a  -  sinh2a*  sinh2p  cos(nj  IN)  , 
for  j  =  1,2,  ....  2N ,  and  a*  satisfies 


(46b) 


sinh2a  sinh2a  *  =  1 . 


(46c) 


Phase  transition  occurs  at  the  critical  temperature  Tc  which  satisfies  the 
equation 


sinh 


(47) 


At  the  critical  temperature  7 ^  =  0.  Following  the  derivation  of  eq.  (46)  in  [63] 
we  can  easily  see  that  yr>0,  for  every  r  =  1,2, .  . . , ZN-l  and  that 
72tv  >0,  if  T  <TC ,  whereas,  y^cO,  ifT>Tc.  Therefore,  eq.  (46),  together  with  the 
identity 

6  We  use  capital  letters  for  x,y,z,  co  to  denote  that  they  take  values  in  (-1,  +1  }.  Small 
letters  will  denote  variables  x,y,z,t u  which  take  values  in  {  0,  +  1 }. 
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i- 1 


cosh/8  =  2l  1  {”[ 


f-0 


cosh8-cos 


(2s  + 1  )k 
21 


and  the  knowledge  of  the  signs  of  y/s,  can  be  used  to  calculate  the  partition 
function  of  the  Ising  model.  In  the  case  of  large  lattices,  and  for  temperatures 
T  *TC ,  the  following  approximation  can  be  used  [61] 

1 


MN 


1  nZMN  =  —  In  ( 4  cosh2a  cosh2(J }  +  B  (a,  ^ )  , 


(48a) 


where 


1  r2lt  r2* 

B  (a,  (3 )  =  — j  Jo  Jo  ln<t>(tanh2a,  tanh2B;0,  t>)d.Qd.t>  , 


(48b) 

(48c) 


<t>(x ,y ; 0, 0)  =  1  -  x  (1-y2)2  cos9  -  y  (1-x2)2  cos<l)  , 

with  x  =  tanh2a  and  y=tanh2(3.  Equation  (48b)  can  be  approximated  by  its 
Riemann  sum;  therefore, 


B  (a,  8 )  =  — ^rr  £  £  In  4>(tanh2a,  tanh2p ;  0;* ,  0*  )  , 

j=\  *=l 


(49a) 


where 


0,-j,  = 


2j  -2  2k -l 


M 


MN 


n  , 


(49b) 


and 


<t>*  = 


2k  -1 


N 


(49c) 


Equations  (48)  and  (49)  are  quite  accurate,  even  in  the  cases  of  small  lattices; 
therefore,  we  shall  use  them  when  M ,  N  >  20.  When  M ,  N  <20,  eq.  (46)  should 
be  used. 

Equations  (48)  and  (49)  can  also  be  used  to  compute  the  first-  and 
second-order  derivatives  of  the  partition  function  with  respect  to  T.  Indeed,  it 
is  easy  to  show  that, 


(50a) 


Emn(T) =  T 


2  1  31  nZ\fN 

mn  ar 


1  ain ZS!N 


1  a2lnZ.v 


whert; 


MN 

?T  + 

M/V  ar2 

1 

ainZMW 

A  1 

c/inZ^'v 

B  1  31nZWN 

MN 

ar 

T2  MN 

3a 

r2  MiV  3p 

1 

32lnZ,w 

1  f  2 A 

ainZ'rv 

2  B  31  nZMN 

MiV 

ar2 

miv  [  r3 

3a 

+  T3  3p 

A5 

32ln Zmn 

2 AB  31n ZMN  B2 

1  , 

rp4 

3a2 

T 4  3a3p  r4 

(50b) 


(50d) 


and,  for  T  *TC, 


1 

d\nZMN 

MN 

da 

1 

31  nZMN 

MN 

3P 

1 

32ln ZMN 

MN 

3a2 

1 

32lnZ.w 

MiV 

3p2 

1 

32ln ZyN 

MiV 

3a3p 

.  in  1  "  £  0«Cxfy;e;*,<>*) 

=  tanh?a  +  V  V  — - — — ~ - —  , 

2 MN  jtl  Atl  <t>(x ,  y ;  9;* ,  0* ) 

^  2MN  Jj  ^  <I>(*  ,y ;  9y* ,  <t»* ) 

-  2(1  -  tanh^a)  >  -Lg.  [g  i 

2MiV  jTuti  <J>* ) 


_  £  £  ^Oe.y;^*,^) 

y=i*=i  ^(x.y ;0/*,<t>*) 

l  £  £  <&!»(» ,y;8j*,»t) 
WiV  ^,y;0;*,o*) 

jtt *r-  <^(x,y;0,*,0*) 


=  2(1-  tanh22(J)  + 


(51a) 

(51b) 


(51d) 


1  *  "  fl(*,y  ;9;*,'D*) 

^  ^  /W*.  A  1 


_  £  £  ^  ,y ;  ,  o* )  <&cl*  ,y ;  8/*  >  a* ) 

h  * = i  ^(x  ,y  ;0j* ,  c>* ) 


(51e) 
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In  eqs.  (51)  Oa,  Cw,  <t>&,  <&»  and  <t>o8  are  the  first-  and  second-order  partial 


derivatives  of  4>  with  respect  to  a  and  {5,  given  by 

_i_  _i_ 

4>a(x,y;Q,<|))  =  -2(l-*2)(l-y2)2  COS0  +  2 xy  (1-x2)2  cos<j>  ,  (52a) 

j_  i_ 

<D(i(x,y  ;6,<>)  = -2(l-y2)(l-x2)2  cos0  +  2xy  (1-y2)2  cos0  ,  (52b) 

(&aa(x,y  ;0,O)  =  8x  (l-*2)(l-y2)2  cosi)  +  4(l-2x2)y  (I-*2)2  cosb  ,  (52c) 

i_  _i_ 

<t>wiGc,y  ;9,0)  =  8y  (l-y2)(l-r2)2  cos0  +  4(l-2y2)*  (1-y2)2  cose  ,  (52d) 

and 

‘t’ad^.y ; 0, 4>)  =  4y  (l-jc2)(l-y2)2  cos0  +  4x  (1-x2)2  (l-y2)cos0  .  (52e) 


Equations  (48)-(52)  are  used  for  the  analytical  computation  of  the  parti¬ 
tion  function  and  its  derivatives.  To  obtain  the  Monte-Carlo  estimates,  we  have 
to  modify  our  program,  since  the  random  variable  assumes  values 

- 1,  + 1  instead  of  values  0,  + 1,  and  since  we  have  assumed  toroidal  boundary 
conditions,  instead  of  free  boundary  conditions.  As  the  lattice  becomes  large 
the  effects  of  the  boundary  conditions  become  negligible  [2]. 7  We  can  now  use 
the  original  Monte-Carlo  simulation  program,  used  in  the  first  class  of  experi¬ 
ments,  by  employing  the  following  transformation 

X  =  2x  -  1  ,  Y=2y-1,  Q  =  2cd-1,  (53) 

which  yields  binary  0-1  random  variables.  Substituting  eq.  (53)  into  eq.  (45)  we 
obtain 

<j(x  ,y  ,z,co)  =  exp  -CA  +  B)-(2A  +  2£)x -2B  y -2Aco  +  43  xy +4A  xcol  . 

7  This  may  not  be  the  case  when  the  temperature  is  close  to  the  critical  temperature,  in 
which  case  a  discrepancy  between  the  analytical  and  computational  result  is  expected. 
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In  our  simulation  experiments  we  have  considered  four  different  Ising 
models  defined  on  various  size  lattices,  up  to  128x128  sites.  These  models  are 
viewed  as  one-parameter  models  in  which  A,  B  are  kept  fixed  and  T  varies. 
Realizations  of  these  GRF’s  on  a  128x128  site  lattice,  at  different  tempera¬ 
tures,  are  depicted  in  Fig.  9.  The  corresponding  values  of  A  and  B ,  together 
with  some  other  useful  information,  are  depicted  in  Table  III.  Our  objective 
now  is  to  compare  the  exact  values  of  the  partition  function  with  the  estimated 
ones,  obtained  by  using  Method  1  and  Method  2.  Figures  10A  and  10B  show 
the  typical  convergence  behavior  of  our  two  methods  for  the  case  of  the  Ising  1 
model  considered  on  a  16x16  site  lattice,  at  a  high  temperature  (weak  cou¬ 
plings)  and  at  a  low  temperature  (strong  couplings).  The  obtained  results  verify 
our  previously  observation  that  Method  2  is  more  accurate  at  temperatures 
below  the  critical  temperature  than  Method  1.  In  Figs.  11A-11D  we  compare 
the  estimates  obtained  by  using  Method  1  and  Method  2  with  the  exact  results 
for  the  four  Ising  models  of  Table  III  considered  on  a  32x32  site  lattice.  The 
same  quantities  are  depicted  in  Figs.  12A-12D,  for  a  64x64  site  lattice,  and  in 
Figs.  13A-13D  for  a  128x128  site  lattice.  In  all  cases  K  -  Kr  =  K2  =  10,000 MN . 
The  obtained  results  indicate  that  Method  2  provides  really  good  estimates  of 
the  partition  function,  whereas,  Method  1  fails  to  do  so  at  low  temperatures.  In 
Figs.  14A-14D  we  compare  E32i32(T ;  Ku  K2)  with  E32t32(T),  and  C32,32{T\ Ku  K2) 
with  C 32,32(T),  in  order  to  check  on  the  accuracy  of  Algorithm  III.  We  have  used 
Kx  =  K2  =  322x  10,000  iterations.  Finally,  we  consider  the  Ising  model  as  a  two- 
parameter  model,  with  respect  to  the  parameters  (a  =  A  IT,  |3  =  B/T),  on  a  16  x  16 
site  lattice.  We  took  0.2  <  a,  (3  <  0.8,  which  corresponds  to  models  with  both 
strong  and  weak  couplings,  as  it  is  indicated  in  Fig.  15.  The  partition  function 
of  the  Ising  model  has  been  calculated  by  using  Method  1  (see  Fig.  16B)  and 
Method  2  (see  Fig.  16C).  The  exact  partition  function  is  depicted  in  Fig.  16A. 
The  resulting  computational  error  is  shown  in  Figs.  17A  and  17B.  These 
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results  clearly  illustrate  the  superiority  of  Method  2. 

To  conclude,  we  would  like  to  point  out  that,  although  Method  2  is  usually 
superior  to  Method  i  at  virtually  all  temperatures,  we  prefer  to  use  Method  1 
at  high  temperatures  (i.e.,  for  T  >TC)  because  of  its  reasonable  accuracy  and 
computational  efficiency  (recall  that  the  implementation  of  Method  2  requires 
the  computation  of  the  first-order  derivatives  of  the  partition  function,  with 
respect  to  cr(x  ,y,z,c o),  which  is  a  quite  expensive  computation). 

VII.  CONCLUSIONS 

We  have  presented  a  new  technique  for  the  estimation  of  the  partition 
function  of  a  general  GRF  which  is  based  on  approximating  a  general  GRF  by 
a  MC-GRF.  We  have  discussed  an  optimal  choice  for  the  MC-GRF  in  terms  of 
achieving  error  variance  reduction,  and  we  have  adopted  two  reasonable  subop- 
timal  choices  that  have  the  advantage  of  simplicity  and  computability.  These 
choices  contain  substantial  information  about  the  initial  GRF,  and  allow  us  to 
achieve  significant  error  variance  reduction.  The  second  choice  is  optimal  in 
terms  of  minimizing  an  entropy  distance  from  the  given  Gibbs  distribution.  We 
have  also  considered  different  Monte-Carlo  algorithms  and  concluded  that  the 
choice  of  the  Gibbs  Sampler  with  lexicographic  site  updating  was  the  most 
appropriate.  Our  proposed  techniques  result  in  unbiased  and  consistent  esti¬ 
mates  of  the  partition  function.  This  is  a  major  improvement  over  existing 
Metropolis-like  simulation  algorithms,  which  sample  directly  from  a  distribu¬ 
tion  that  asymptotically  approaches  the  Gibbs  distribution,  and  are  capable  of 
estimating  only  the  derivatives  of  the  logarithm  of  the  partition  function,  and 
not  the  logarithm  of  the  partition  function  itself. 

Our  two  choices  of  MC-GRF’s  have  resulted  in  two  different  methods  for 
the  calculation  of  the  partition  function.  We  have  carried  out  many  simulation 
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experiments  in  order  to  test  the  reliability,  accuracy  and  the  relative  merits  of 
the  two  approaches.  The  obtained  results  have  been  compared  to  analytical 
ones,  wherever  possible.  Our  experiments  gave  remarkably  accurate  results, 
demonstrating  the  fact  that  we  can  obtain  good  estimates  within  reasonable 
computational  time,  for  a  wide  variety  of  GRF  models.  We  have  also  examined 
the  behavior  of  the  methods  as  the  lattice  size  increases  and  as  the  tempera¬ 
ture  approaches  the  critical  temperature,  and  conjectured  that  Method  2  is 
more  reliable  than  Method  1  around  and  below  the  critical  region  of  the  GRF. 
Nevertheless,  we  prefer  to  use  Method  1  at  high  temperatures.  We  believe  that 
our  methods  can  be  effectively  used  for  the  optimal  parameter  estimation  of  a 
general  GRF.  We  would  finally  like  to  notice  that,  the  use  of  computationally 
efficient  algorithms  can  boost  the  performance  of  our  methods,  and  yield  even 
more  reliable  estimates. 
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Figure  1:  Realizations  of  five  128x128  site  GRF’s  at  different  temperatures.  Each 
row  depicts  a  GRF  model  (top  to  bottom:  GRF1,  GRF2.  GRF3.  GRF4. 
GRF5)  at  five  different  temperatures  'left  to  tight:  7\,  Tu  7ti,  T  ■,  T 
See  also  Tables  1,  II. 
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Table  I:  Temperatures  used  in  the  GRF 
realizations  depicted  in  Fig.l. 


Reference 

name 

To 

Temperatures 

T  i  T2  T3 

T 4 

GRFl 

2.00 

1.30 

1.00 

0.75 

0.50 

GRF2 

2.00 

1.70 

1.40 

1.20 

1.00 

GRF2 

3.00 

2.00 

1.30 

1.00 

0.50 

GRF4 

3.00 

2.00 

1.30 

1.00 

0.80 

GRF5 

3.00 

2.30 

2.00 

1.30 

1.00 

— 

Table  II:  LTF’s  of  the  five  GRF’s  depicted  in  Fig.l. 


LTFn 

GRFl 

GRF2 

GRF3 

GRF4 

GRF5  i 

T  x  lna(Q,Q,0,Q) 

0.00000 

0.00000 

0.0QQQQ 

0.00000 

0.00000 

T  x  In  a(0, 0,0,1) 

-0.50000 

0.00000 

0.00000 

-0.40000 

0.00000 

T  x  lna(0,0,l,0) 

-0.50000 

0.00000 

0.00000 

0.03000 

0.00000 

T  x  In  a(0, 0,1,1) 

-1.00000 

0.00000 

0.00000 

-0.72000 

0.00000 

T  x  In  a(0, 1,0,0) 

1.00000 

0.00000 

0.00000 

-0.91800 

0.00000 

T  x  In  a(0, 1,0,1) 

0.51500 

-1.25000 

0.59500 

-2.04300 

5.00000 

Tx  In a(0, 1,1,0) 

0.50000 

0.00000 

0.00000 

-0.55000 

0.00000 

T  x  In  a(0, 1,1,1) 

0.01500 

-1.84960 

0.59500 

-2.85500 

5.20000 

Tx  In  oil,  0,0,0) 

-0.26000 

0.89540 

0.00000 

-3.20000 

0.20000 

Tx  In o(l, 0,0,1) 

-2.76000 

-0.94740 

-1.25000 

-0.60000 

-1.80000 

Tx  In o(l, 0,1,0) 

-0.63000 

1.22990 

0.59500 

-4.99000 

-4.30000 

Tx  In o(l, 0,1,1) 

-3.13000 

-0.68070 

-0.65500 

-0.24000 

-6.70000 

Tx  In a(l,l, 0,0) 

2.84000 

-0.86980 

-1.25000 

-0.64260 
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Figure  2A 

GRF1.  4x4  sites 
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GRF2,  4x4  sites 
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Figure  2C 

GRF3,  4x4  sites 
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Figure  2E 
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Figures  2A-2E:  Comparison  of  the  exact 
and  estimated  partition  function,  negative 
internal  energy  and  specific  heat  for  the 
five  GRF  models  of  Table  II  (K  = 
Kx=K2  =  1600 ). 
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Figure  4A 

GRF1,  16  x  16  sites 
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GRF2,  16  x  16  sites 
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Figure  4C 

GRF3,  16  x  16  sites 
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Figure  4E 

GRF5,  16  x  16  sites 
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Figures  4A-4E:  Monte-Carlo  estimation 
of  the  partition  function  of  the  five  GRF 
models  of  Table  II,  using  three  different 
sampling  schemes  (K  =  162x  io4). 
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GRF3,  64  x  64  aitea 
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Figure  6D 
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Figure  6E 

GRF5,  64  x  84  aitea 
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Figures  6A-6E:  Monte-Carlo  estimation 
of  the  partition  function  of  the  five  GRF 
models  of  Table  n,  using  three  different 
sampling  schemes  (K  =  642x  10*). 
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Figure  7A 

GRF1.  128  x  128  sites 


Figure  7B 

GRF2,  128  x  128  sites 
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Figure  7C 
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Figure  7E 

GKF5,  128  x  128  sites 
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Figures  7A-7E:  Monte-Carlo  estimation 
of  the  partition  function  of  the  five  GRF 
models  of  Table  II,  using  three  different 
sampling  schemes  (K  =  1282x  104). 
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Figure  9:  Realizations  of  128x128  site  Ising  models.  Each  row  depicts  an  Ising 
model  (from  top  to  bottom:  Ising  1,  Ising  2,  Ising  3,  Ising  4)  at  five 
different  temperatures  (from  left  to  right:  T0,  Tu  T2,  T2,  Tj.  See  Table 
III  for  details. 
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Table  HI:  Interaction  parameters  04,5)  of  the  Ising  models 
of  Fig.  9  for  different  temperatures  T, . 


Reference 

name 

Interactions 
(A  ,5) 

Crit.  temp. 

T 

Temperatures  of  Fig.  9 

T0  Tt  t2  r3  t4 

Ising  1 

(1.0, 1.0) 

2.26918 

4.00 

2.27 

IPSM 

Ising  2 

(2.0, 1.0) 

3.28204 

4.00 

3.28 

3.00 

Ising  3 

(3.0, 1.0) 

4.15617 

6.00 

4.15 

3.50 

Ising  4 

(-1.0, 1.0) 

2.26918 

4.00 

2.50 

1.50 

Figure  10A 

ISING  1,  16  x  16  aitea,  T  =  3.0 


Iterationa 

Figure  10B 

ISING  1,  16  x  16  aitea,  T  =  1.5 


Figures  10 A,  10B:  Convergence  behavior  of  the  Monte-Carlo  algorithm 
for  the  estimation  of  the  partition  function  of  the  16x16  site  Ising  1 
model  at  T  =  3.0  (Figure  10A,  T  > Tc ),  or  at  T  =  1.5  (Figure  10B,  T  <TC). 
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Figures  11A-11D:  Monte-Carlo  estimation  of  the  partition  function  of  the  four  Ising 
models  of  Table  III,  considered  on  a  32  x  32  site  lattice,  by  using  two  different  sampling 
schemes  ( K =322xl04).  Comparison  with  the  exact  partition  function. 

& - — K-lnZ32  32  -  Exact  result. 

322 


^lnZJ23Zp.(K)  -  Method  1. 
^7  In  *»*,-«*)  -Method  2. 


UK)  -  Method  2. 


page  67 


Figure  16B 


Figure  16C 


Figures  16A-16C:  Comparison  of  the  exact  and  the  Monte-Carlo  estimated  partition 
Junction  of  a  16x16  site  Ising  model  with  parameters  0.2<a<0.8  and  0.2<3<0.8 
(K  =  162x  104). 
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Figures  17A,  17B:  Error  (%)  of  the  Monte-Carlo  estimate  of  the  partition  function  for 
the  Ising  model  considered  in  Figs.  16A-C. 
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